% Obtains listing-level surplusses given commissions and valuation
% parameters. Returns matrices with in the first dimension the seller's
% valuation, and in the second dimension the number of bidders. Corresponds
% to equations 4-7 in the paper. 
% Run this first before getting to v0star, evaluate this once and then use
% it in subsequent bidder and seller surplus calculations.
% ---------------------------------------------- %
function[pibmat,pib0vec,pdfmat,pismat,pis0mat,u0vec,rvec,pipmat,pip0vec,pricemat,highvalmat,reservemat,saleprobmat,price0vec,highval0vec,volumemat,volume0vec]...
    = piBS_MC(cB,cS,thetab,thetas,PAR,biddervaluethreshold)

    if PAR.generalized==0
        db  = makedist(PAR.dist,thetab(1),thetab(2));
        CDFb = @(x) cdf(db,x);
        PDFb = @(x) pdf(db,x);
        ds  = makedist(PAR.dist,thetas(1),thetas(2));
        CDFs = @(x) cdf(ds,x);
        PDFs = @(x) pdf(ds,x);
        
        % simulate seller tastes
        u0vec_c = unique(sort(icdf(PAR.dist,PAR.v0probs,thetas(1),thetas(2))))';
        u0vec   = unique(sort(icdf(PAR.dist,PAR.v0probsfine,thetas(1),thetas(2))))';

         % extreme values
        PAR.maxv = fminbnd(@(x) (CDFb(x)-0.999)^2,thetab(1)-3,thetab(1)+3);
        PAR.minv = fminbnd(@(x) (CDFb(x)-0.001)^2,thetab(1)-3,thetab(1)+3);
      

    else %Generalized gaussian
        
        % for bidder selection exercise
        if nargin == 6
            CDFb = @(x) (ggd_cdf(x,thetab)-ggd_cdf(biddervaluethreshold,thetab))./(1-ggd_cdf(biddervaluethreshold,thetab));
            PDFb = @(x) ggd_pdf(x,thetab)./(1-ggd_cdf(biddervaluethreshold,thetab));
        else
            CDFb = @(x) ggd_cdf(x,thetab).*ones(size(x));
            PDFb = @(x) ggd_pdf(x,thetab).*ones(size(x));
        end
        
        % extreme values depending on support of FV
        PAR.maxv = fminbnd(@(x) (CDFb(x)-0.999)^2,thetab(1)-3,thetab(1)+3);
        PAR.minv = fminbnd(@(x) (CDFb(x)-0.001)^2,thetab(1)-3,thetab(1)+3);
        location=thetab(1);
        scale=thetab(2);
        shape=thetab(3);
        if shape > 0
             maxv = location+(scale/shape)-0.001;
             PAR.maxv=min(PAR.maxv,maxv);
        else %shape<0
             minv = location+(scale/shape)+0.001; 
             PAR.minv=max(PAR.minv,minv);
        end
        CDFs = @(x) ggd_cdf(x,thetas).*ones(size(x));
        PDFs = @(x) ggd_pdf(x,thetas).*ones(size(x));
        
        % Support
        if thetas(3)>0
            lb = -3*thetas(1);
            ub = min(10,thetas(1)+thetas(2)/thetas(3));
        else
            lb = max(-10,thetas(1)+thetas(2)/thetas(3));
            ub = 3*thetas(1);
        end

        % simulate seller tastes
        u0vec_c = unique(sort(arrayfun(@(i) fminbnd(@(x) (CDFs(x)-PAR.v0probs(i))^2,lb,ub),1:length(PAR.v0probs))));
        u0vec   = unique(sort(arrayfun(@(i) fminbnd(@(x) (CDFs(x)-PAR.v0probsfine(i))^2,lb,ub),1:length(PAR.v0probsfine))));

    end
     
    % Solve for optimal r on a grid
    for i = 1:length(u0vec_c)
        rvec(i) = optreserve(u0vec_c(i),cS,cB,PDFb,CDFb,PAR);
    end    
       
    % Prep u0, r, n matrices
    u0mat    = repmat(u0vec_c,PAR.maxn,1); %v0mat(i,j) is n=1, v0 (or v0bar)=j
    rmat     = repmat(rvec,PAR.maxn,1);
    nmat     = repmat(linspace(1,PAR.maxn,PAR.maxn)',1,length(u0vec_c));
                      
    % Evaluate listing-level (pi) surplus on matrices
    [HVsmall,SHVsmall]   = sampleHSHvalues(PDFb,CDFb,PAR);          % draw values of HV and SHV|HV <- this is v slow but accurate
    HV       = zeros(PAR.maxn,length(PAR.MC));
    SHV      = HV;
    for n = 1:PAR.maxn
        HV(n,:)       = interp1(PAR.probsMC,HVsmall(n,:),PAR.MC,'linear');
        SHV(n,:)      = interp1(PAR.probsMC,SHVsmall(n,:),PAR.MC,'linear');
    end
    pibmata      = rmat;
    pismat_c     = rmat;
    pipmata      = rmat;
    pricemata    = rmat;
    highvalmata  = rmat;
    reservemata  = rmat;
    saleprobmata = rmat;
    volumemata   = rmat;
    
    for i = 1:length(rmat(rmat>-Inf))
        [pibmata(i),pismat_c(i),pipmata(i),pricemata(i),highvalmata(i),reservemata(i),saleprobmata(i),volumemata(i)] = ...   
        EBSsurplus(nmat(i),u0mat(i),rmat(i),cS,cB,HV,SHV,length(PAR.MC),CDFb,PAR);
    end
    
    pdfmat_c   = repmat((arrayfun(@(u) PDFs(u), u0vec_c)),PAR.maxn,1);
    for i = 1:PAR.maxn
        pdfmat_c(i,:)   = pdfmat_c(i,:)/sum(pdfmat_c(i,:));
    end
    pibmat_c = pibmata.*pdfmat_c;  %Buyers need to take expectation over v0 given v0*
    pipmat_c = pipmata.*pdfmat_c;  %Platform needs to take expectation over v0 given v0* - 
    
    % Also evaluate winning bidder surplus and components for welfare analysis
    pricemat_c     = pricemata.*pdfmat_c;
    highvalmat_c   = highvalmata.*pdfmat_c;
    reservemat_c   = reservemata.*pdfmat_c; %for lots that sell
    saleprobmat_c  = saleprobmata.*pdfmat_c;
    volumemat_c    = volumemata.*pdfmat_c;
    
    % When sellers set 0 reserve price
    pib0mata     = nmat;
    pis0mat_c    = pib0mata;
    pip0mata     = pib0mata;
    price0mata   = pip0mata;
    highval0mata = price0mata;
    volume0mata  = price0mata;
    
    for i=1:length(nmat(nmat>-Inf))
        [pib0mata(i),pis0mat_c(i),pip0mata(i),price0mata(i),highval0mata(i),volume0mata(i)]...
        = EBSsurplus(nmat(i),u0mat(i),0,cS,cB,HV,SHV,length(PAR.MC),CDFb,PAR);
    end
 
    pib0vec      = pib0mata(:,1);
    pip0vec      = pip0mata(:,1);
    price0vec    = price0mata(:,1);
    highval0vec  = highval0mata(:,1);
    volume0vec   = volume0mata(:,1);  
    
    % Interpolate on a finer grid
    %pibmat,pismat,pis0mat,pipmat,pdfmat
    pibmat = zeros(length(pibmat_c(:,1)),length(u0vec));
    pismat = pibmat;
    pdfmat = pibmat;
    pipmat = pibmat;
    winbiddermat = pibmat;
    pricemat = pibmat;
    highvalmat = pibmat;
    reservemat = pibmat;
    saleprobmat = pibmat;
    pis0mat = pibmat;
    volumemat = pibmat;         
    for n = 1:PAR.maxn
        pibmat(n,:)  = interp1(u0vec_c,pibmat_c(n,:),u0vec,'linear');
        pipmat(n,:)  = interp1(u0vec_c,pipmat_c(n,:),u0vec,'linear');
        pdfmat(n,:)  = interp1(u0vec_c,pdfmat_c(n,:),u0vec,'linear');
        pismat(n,:)  = interp1(u0vec_c,pismat_c(n,:),u0vec,'linear');
        pis0mat(n,:) = interp1(u0vec_c,pis0mat_c(n,:),u0vec,'linear'); 
        winbiddermat(n,:) = pibmat(n,:)*n;
        pricemat(n,:) = interp1(u0vec_c,pricemat_c(n,:),u0vec,'linear');      
        highvalmat(n,:) =interp1(u0vec_c,highvalmat_c(n,:),u0vec,'linear');
        saleprobmat(n,:) = interp1(u0vec_c,saleprobmat_c(n,:),u0vec,'linear');
        reservemat(n,:) = interp1(u0vec_c,reservemat_c(n,:),u0vec,'linear');   
        volumemat(n,:) = interp1(u0vec_c,volumemat_c(n,:),u0vec,'linear');             
    end        
end
